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Functional magnetic resonance imaging (fMRI) techniques have con- 
tributed significantly to our understanding of brain function. Current 
methods are based on the analysis of gradual and continuous changes 
in the brain blood oxygenated level dependent (BOLD) signal. De- 
parting from that approach, recent work has shown that equivalent 
results can be obtained by inspecting only the relatively large ampli- 
tude BOLD signal peaks, suggesting that relevant information can 
be condensed in discrete events. This idea is further explored here 
to demonstrate how brain dynamics at resting state can be captured 
just by the timing and location of such events, i.e., in terms of a 
spatiotemporal point process. As a proof of principle, we show that 
the resting state networks (RSN) maps can be extracted from such 
point processes. Furthermore, the analysis uncovers avalanches of 
activity which are ruled by the same dynamical and statistical proper- 
ties described previously for neuronal events at smaller scales. Given 
the demonstrated functional relevance of the resting state brain dy- 
namics, its representation as a discrete process might facilitate large 
scale analysis of brain function both in health and disease. 
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Important efforts to understand brain function, both in 
health and disease, are concentrated in the analysis of large- 
scale spatiotemporal patterns of brain activity available from 
fMRI techniques [H El Bl H El E], allowing for instance the 
unravelling of the functional connectivity between all possi- 
ble brain regions, as is done under the Connectome project 
[3 (SI [9] . Novel techniques of analysis are needed because the 
difficulty of managing extremely large data sets is increased by 
new advances in imaging technology continuously improving 
temporal and spatial resolution. In that direction recent work 
has shown that important features of brain functional connec- 
tivity at rest can be computed from the relatively large ampli- 
tude BOLD fluctuations 10 after the signal crosses some am- 
plitude threshold. Pursuing the same general idea, it would 
be desirable to explore complementary ways of data reduc- 
tion. In particular we are interested in a method often used 
to study the structure and properties of attractors of dynam- 
ical systems, which consists in the introduction of a Poincare 
section. By definition, this approach decreases the dimension 
of the phase space, and consequently the size of the data sets, 
facilitating in this way further numerical investigations. In 
general, there exist two possibilities: the first one is to analyze 
the set of points which are the coordinates of the successive 
intersections of the secant Poincare plane by the phase space 
trajectories. The second possibility is to study the series of 
time intervals between the consecutive intersections. The re- 
sulting time intervals constitute a so-called point process ^J , 
a construction useful in many areas of science, including neu- 
roscience. It has been shown that under certain conditions the 
most important statistical features of the dynamical regime 
can be condensed into a point process 12, 13^, 1141 [T5l 1161 flT] . 
The motivation to attempt similar approach in fMRI data 
is strengthened by the observation that, in response to neu- 
ronal activation, the BOLD signal often repeats a stereotypi- 
cal pattern [JOl [THl [HI [20] . This feature of the BOLD response 



dynamics suggests that it should be possible to compress the 
data sets using the temporal marks of some Poincare section 
of the BOLD signal. This is the hypothesis explored here, 
which implies that, in principle, the entire brain resting state 
functional connectivity can be reconstructed solely on the ba- 
sis of the time and location of the BOLD signal threshold 
crossings. Besides its practical importance for fMRI signal 
processing, this approach may provide further clues on the 
dynamic organization of the brain RSN. 

The paper is organized as follows: first the definition of the 
point process is sketched, followed by its validation by repli- 
cating well studied aspects of the fMRI brain resting state. 
The spatiotemporal statistics are then considered, revealing 
novel aspects of the brain dynamics which are scale invariant, 
consistent with that shown for other systems at the critical 
state L2ll|22l[23l[2l]. 
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Fig. 1. (A) Example of the normalized BOLD signal of one voxel and its point pro- 
cess (filled circles) defined by the threshold (dashed line at 1 S.D.) crossings events. 
(B) Average BOLD signal (for one subject) triggered at each threshold crossing as in 
10. 



Results 

The fMRI dataset is reduced to a spatiotemporal point pro- 
cess by first normalizing each BOLD signal by its own S.D., 
and subsequently selecting the time points at which the in- 
creasing signal crosses a given threshold (1 S.D. in this case) 
as it is shown in the example of Figure 1. Notice that the av- 
erage BOLD signal around the extracted points (Figure IB) 
resembles the typical hemodynamic function seen in response 
to an stimulus ^18. .19^ , despite the fact that in this case there 
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are not explicit inputs, since recordings were obtained under 
resting conditions. 

The timing and spatial location of the extracted points 
is the only information needed for further analysis. For the 
parameter used here, from each voxel BOLD time series (240 
samples) only 15 ±3 points are threshold crossings (about one 
point every 40 sec) which corresponds to near 94% reduction 
of the data, (see additional estimations in Supp. Info). 
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Fig. 2. RSN maps constructed with the point process compared with standard 
PICA. (A) PICA spatial maps (left column) and rate of points conditional to activity 
at a given seed (rightmost three columns, each one corresponds to a different seed). 
(Slice z coordinates are -12, 0, 0, 36, 20, 26 for RSN 1 to 6; for coordinates see 
Table SI). Scales for PICA (Zpjca) and conditional rate {Zcr) calculations are 
depicted in the inset. (B) Conditional rate maps constructed using 3, 6 and 12 events 
of the point process at the ANGL seed (averaged from ten subjects. Slice coordinates 
are x=-4,y=-60,z=18) (C) Correlation between RSN5 (the default-mode network, 
DMN) PICA-derived map and the point process-derived conditional rate maps, as a 
function of the number of points used. Arrows denote the examples of panel B. Z 
scores (number of points as degrees of freedom) with the line of 95% confidence are 
plotted in the inset. 



Resting state networks maps derived from the 
point process. Despite the very large data reduction, we 
found that the information content of the few remaining points 
is very high. To demonstrate this we use the point process 
to calculate the location of six well known RSN maps. The 
computation of these maps are compared with those obtained 
from the full BOLD signal using a well established method 
(probabilistic independent component analysis - PICA [6]). 
This is done by calculating in six RSNs the rate of points 
co-occurrence (up to 2 time units later in this case) between 
representative sites ("seeds") and all other brain voxels and 
presented as maps in Figure 2. The seeds locations were se- 
lected according with previous work (see coordinates in Table 
SI). The similarities between our conditional rate maps and 
the respective PICA maps (rightmost three columns and left 
column of Figure 2A respectively) is already obvious to the 
naked eye and confirmed by the correlation plotted in Fig. 
2C. The calculation shows that despite using less than 6% of 
the raw fMRI information, about 5 points (on average) are 
enough to obtain RSN maps that are highly correlated (95% 
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Fig. 4. Typical fluctuations of the clusters. (A) Examples of co-activated clus- 
ters of neighbors voxels (depicted in different colors). (B) Example of the temporal 
evolution of the clusters' number and maximum size (in units of voxels) in one indi- 
vidual. (C) Instantaneous relation between the number of clusters vs. the number 
of active sites showing a positive/negative correlations depending whether activity is 
below/above a critical value (~ 2500 voxels, indicated by the dashed line here and 
in Panel B). (D) The cluster size distribution, follows an inverse power law spanning 
four orders of magnitude. Individual statistics for each of the ten subjects are plotted 
in the inset. 



confidence) with those obtained using PICA of the full BOLD 
signals. 

Structure and dynamics of the active clusters. The 

analysis of the brain large scale spatiotemporal patterns using 
the full BOLD signals is hampered by numerical constraints. 
In contrast, the analysis of brain dynamics with the derived 
point process is simpler by compiling the statistics of clus- 
ters of points, both in space and time. Clusters are groups of 
contiguous voxels with signal above the threshold at a given 
time, identified by a scanning algorithm in each fMRI vol- 
ume. Figure 3A shows examples of clusters (in this case non- 
consecutive in time) depicted with different colors. Typically 
(Fig. 3B top) the number of clusters at any given time varies 
only an order of magnitude around the mean (^ 50). In con- 
trast, the size of the largest active cluster fluctuates widely, 
spanning more than four orders of magnitude. 

The analysis reveals four novel dynamical aspects of the 
cluster variability which hardly could has been uncovered with 
previous methods. 1) At any given time, the number of clus- 
ters and the total activity (i.e., the number of active voxels) 
follows a non-linear relation resembling that of percolation 
[2T [. At a critical level of global activity {^ 2500 voxels, 
dashed horizontal line in Fig. 3B, vertical in Fig. 3C) the 
number of clusters reaches a maximum (~ 100 — 150) as well 
as its variability. 2) The correlation between the number of 
active sites (and index of total activity) and the number of 
cluster reverses above a critical level of activity, a feature al- 
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ready described in other complex systems anytime that some 
increasing density competes with hmited capacity [21 . 3) 
The rate at which the very large clusters (i.e., those above 
the dashed line in 3B) occurs (^ one every 30-50 sec) corre- 
sponds to the low frequency range at which RSN are typically 
detected using PICA [51 [6] . 4) The distribution of sizes (Fig- 
ure 3D) reveals a scale free distribution. These four features 
are known to be generic properties of other complex systems 
EH ^2\ |23 and as such it should be further explored wether 
they are preserved under developmental changes, aging and 
pathological conditions. Notice, for instance, that the non- 
linear correlation (i.e. Fig. 3 B and C) between the number 
of clusters and the size of the largest activated cluster can 
be related (and an objective measure) to the balance between 
excitation and inhibition, a property believed to be funda- 
mental for healthy cortical function. In the same direction, 
the size of the largest cluster and the number of clusters can 
be seen in terms of the amount of integration and segregation 
long advocated by Tononi et al. [25l [26] as the fundamen- 
tal conundrum that the healthy cortex needs to be executing 
at any given time. The important issue to be noted here is 
that the point process approach allows for the first time a 
straightforward quantification of these properties. 
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Fig. 3. Clusters spread throughout the brain as scale-free avalanches. (A) Two 
examples of avalanches, one triggered from the visual cortex (top) and another from 
insular cortex (bottom). (B) Average cluster's fractal dimension D ~ 2.15 dz 0.02 
estimated by the slope of the counts vs. length plot. Inset: derivatives between points 
in the main plot for each subject. (C) Avalanche's size and lifetime distribution func- 
tion computed from about 8000 avalanches in each of 10 subjects. While avalanche 
size follows an inverse power law, their lifetimes density decrease faster. The dashed 
line corresponds to an exponent of 3/2, found previously for neuronal avalanches |27j . 



Activity spread is scale-free. Now we turn to explore 
in detail how clusters evolve in space. This is motivated by the 
usual notion that cognition and behavior underlying mecha- 
nism involves the sequential activation of clusters of coher- 
ent events through regions of the cortex. In that sense, we 
track over time an activated cluster which can appear, grow 
to achieve a maximum size and then disappear (or translate 



or divide into sub-clusters) as shown in the example of Fig. 
4A. the data reduction allow us to characterize easily this pro- 
cess, and in particular study two properties. For each cluster, 
first we measured a static space filling property, the average 
fractal dimension D. This is shown in Fig. 4B which illus- 
trate that D ^ 2.15 =b 0.02. Second we look at the dynamics 
of the cluster propagation, which we found happens in bursts. 
The statistics in Fig. 4C shows that avalanche's could last up 
to 30 sec. with sizes up to 10^ and have no preferred scale, 
being fitted to an inverse power law with the same 3/2 expo- 
nent described previously in smaller scales [27, 28, 23 . When 
we look at the spatial coverage of typical avalanches at their 
maximum size, we found that is not arbitrary, because they 
remain confined to the RSN' contiguous regions at which they 
started. For instance, in the RSN5 (i.e., DMN) an avalanche 
could start at the Medial Prefrontal Cortex (mPFC) would 
cover the entire mPFC, but will not extend to the Angular 
Gyrus which is consistent with the notion of RSN being inde- 
pendent components (see Supp. Info.). 



Discussion 

As far as we know, this is the first attempt to describe large 
scale brain dynamics as a point process. The only previous 
report we are aware of 29 dealt with the reverse process: 
how to model the continuous fMRI signal starting from a spa- 
tiotemporal point process. 

Concerning the reason behind the surprisingly successful 
data reduction achieved, we can only conjecture that the up- 
ward going BOLD signals are nonlinear events where the times 
of crossings preserve the most relevant information. This 
could be in line with recent findings of all-or-none "coher- 
ence potentials", exhibiting avalanches with identical scaling 
behavior [3^ in monkey cortex. In a similar fashion, the re- 
duction of the fMRI BOLD signal to discrete events not only 
allows the identification of well-described resting state net- 
works but also it is shown to organize in neuronal avalanches. 

The fact that avalanches with the same power law are 
observed at a large scale is consistent with the hypothesis 
that brain dynamics operates at a critical point of a second 
order phase transition [23], thus rendering the system scale 
invariant and allowing its description at all scales by the same 
set of statistical relations. The spatio-temporal scale invari- 
ance shown in the findings of this work and previous electro- 
physiological experiments on neuronal avalanches is further 
supported by the observation that the correlation function of 
fMRI BOLD signals exhibits fractal properties [24 and the 
correlation length of activity measured with fMRI diverges as 
predicted by the theory of phase transitions [31 . 

The observation that large scale brain dynamics can be 
traced as discrete events (scale free avalanches of activity) 
raises the question of the physiologically relevance encoded 
in the timing of these events. For instance, although rare, 
avalanches in the tail of the power law distribution emerge 
from a local origin and propagate as far as the length of the 
entire cortex, suggesting a role in the binding processes of far 
apart cortical regions. Additionally, the nonlinear relation be- 
tween activated cortical tissue and number of clusters exhibits 
an optimal point, in which the highest levels of brain activity 
are segregated into the maximum number of spatially isolated 
activations. It is interesting to investigate whether total or 
partial disruption of these large events, as well as alterations 
in the balance between activation and segregation into clus- 
ters could be correlated with pathological conditions and with 
the level of awareness of the subject. 

The present results are consistent with observations at 
faster time scales, in which the scalp EEC is reduced to a cer- 
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tain number of stereotypical topographical maps (i.e., EEG 
microstates) j32] and the observation of non-stationarities 
with define discrete segments of electrical activity ^J . Both 
descriptions of electrical scalp activity have also been shown 
to exhibit properties consistent with critical dynamics J34ll35] . 
Further multimodal imaging studies could link these observa- 
tions together in the context of discrete avalanches of neural 
activity propagating through the cortex and determine their 
functional relevance for health and disease. 

In summary, the results show that the location and timing 
of the largest BOLD fluctuations define a spatial point process 
containing substantial information of the underlying brain dy- 
namics. Despite the very large data reduction (> 94%), the 
approach was validated by the favorable comparison of the 
conditional rate maps of avalanching activity with those con- 
structed with the fuU fMRI BOLD signals using PICA. In ad- 
dition to uncover new dynamical properties for the activated 
clusters the method exposed scale- invariant features conjec- 
tured in the past 23 which are identical to those seen at 
smaller scales 27, 28, 23 . Beyond its potential value for fMRI 
signal processing, the ability of the present approach to cap- 
ture relevant brain dynamics underscoring nonlinear aspects 
of the BOLD signal deserve further exploration. 

Material and Methods 

fMRI data acquisition and preprocessing. Data was ob- 
tained, after informed consent, from ten right-handed healthy 
volunteers (9 female, 1 male; mean age=49, S.D.=12) scanned 



each 2.5 sec. during 10 minutes, requested to keep their eyes 
closed and to avoid falling asleep (see Supp. Info for scan- 
ning parameters). The study was approved by the Clinical 
Research Ethics Committee of the University of the Balearic 
Islands, (Palma de Mallorca, Spain). The Melodic package 
was used for the calculation of ICA of the RSN [5] in (Fig. 2) 
as well as for de-noising motion artifacts. 

Cluster and avalanche analysis. Spatial clusters of ac- 
tivated voxels were identified using an algorithm implemented 
in MATLAB, based on the detection of connected compo- 
nents in a CO- activated first neighbors graph (see Supp. Info 
for a description of the algorithm). Clusters' fractal dimen- 
sion was calculated using a standard box-counting algorithm. 
Avalanches were defined (similarly to that done in sandpile 
models, and others [21] [22]) as starting with the isolated ac- 
tivation (i.e., not by any of its neighbors) of a previously in- 
active voxel (or group of voxels), continuing while at least 
one contiguous voxel is active in the next time step and oth- 
erwise ends. The avalanche tracking algorithm implemented 
in this work uses as a criteria for avalanche membership non 
empty intersection with a previously identified cluster of the 
avalanche at a previous time. This is able to resolve shrinking 
and expanding of clusters, translation and division, whenever 
there is spatial overlap at subsequent times. For a technical 
description of the algorithm, see Supp. Info. 
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